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Stability of solitons in parity-time (PT)-symmetric periodic potentials (optical lattices) is analyzed 
in both one- and two-dimensional systems. First we show analytically that when the strength of 
the gain- loss component in the PT lattice rises above a certain threshold (phase-transition point), 
an infinite number of linear Bloch bands turn complex simultaneously. Second, we show that while 
stable families of solitons can exist in PT lattices, increasing the gain-loss component has an overall 
destabilizing effect on soliton propagation. Specifically, when the gain-loss component increases, 
the parameter range of stable solitons shrinks as new regions of instability appear. Thirdly, we 
investigate the nonlinear evolution of unstable PT solitons under perturbations, and show that the 
energy of perturbed solitons can grow unbounded even though the PT lattice is below the phase 
transition point. 
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I. INTRODUCTION 

Recent interest in study of parity-time (PT) symmetric 
optical potentials has its roots in quantum mechanics. In 
quantum mechanics, in order for the energy levels to be 
real and the theory to be probability conserving, it is usu- 
ally assumed that the Hamiltonian (Schrodinger) opera- 
tor be Hermitian. However, in the past decade there has 
seen considerable attention [H-[3l in a weaker version of 
the Hermiticity axiom which requires that the Hamilto- 
nian instead only exhibit space-time reflection symmetry 
(PT symmetry). While there has been much theoreti- 
cal success in developing a non-Hermitian quantum field 
theory, the phenomena unique to this class of pseudo- 
Hcrmitian systems have not yet been observed experi- 
mentally. 

The same Schrodinger equation from quantum me- 
chanics applies also to optics. Motivated by this con- 
nection, optical systems which have PT-symmetric po- 
tentials have been formulated A PT-symmetric op- 
tical potential V(x) is realizable by the careful distri- 
bution of gain and loss in the media so that it satisfies 
the PT symmetry V(x) = V*(— x), where x is the spa- 
tial coordinate and '*' stands for complex conjugation. 
That is, the refractive-index profile of the media is even 
and gain-loss profile is odd. Such optical PT media have 
been created experimentally @, 0] ■ These linear PT me- 
dia undergo phase transition as the gain-loss component 
crosses a certain threshold 0, 041] ■ Below this threshold, 
all eigenvalues of the PT potential are real; but above 
this threshold, complex eigenvalues appear, hence the in- 
tensity of a light beam grows exponentially during linear 
propagation. The nature of this phase transition (es- 
pecially for periodic PT potentials) has not been fully 
understood yet. 

These phenomena may also be studied in a nonlinear 
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context by considering the existence of localized modes 
called solitons [8|, |9[ . When a system contains gain and 
loss, solitons generally exist only at special values of 
the propagation constant [l(J. However, since PT po- 
tentials can admit all-real linear spectra, solitons could 
exist at continuous ranges of the propagation constant, 
which is quite remarkable. So far, soliton families in PT- 
symmetric periodic potentials with defects and in PT- 
symmetric nonlinear potentials have been investigated 
[Ill4l4j . But stability properties of these PT solitons (es- 
pecially in periodic PT potentials) have not been care- 
fully examined. 

In this paper, we investigate linear phase transition 
and stability of (nonlinear) solitons in PT-symmetric pe- 
riodic potentials (optical lattices) in both one and two 
spatial dimensions. Our mathematical model is the non- 
linear Schrodinger (NLS) equation with a PT lattice po- 
tential, 

iU, + U xx + U yy + V(x, y)U + a\U\ 2 U = 0, (1.1) 

where a — ±1 denotes the focusing and defocusing non- 
linearity, and the potential V{x,y) is periodic in x and 
y and satisfies the PT symmetry V(x,y) — V*(— x, —y). 
For simplicity, we take this PT lattice potential to be 

V(x) = V [cos 2 {x) + iW sm{2x)] (1.2) 

in one dimension (ID) and 

V(x, y) = V {cos 2 (x) + cos 2 (y) + iW [sia(2a;) + sin(2y)]} 

(1.3) 

in two dimensions. Here Vo (> 0) is the depth of the real 
component of the potential, Wq is the relative magnitude 
of the imaginary component, and the period of this PT 
lattice is tt. For this system, we first show analytically 
that when the strength of the gain-loss component (the 
imaginary part of V(x, y)) in the PT lattice rises above 
a certain threshold (phase-transition point), an infinite 
number of linear Bloch bands turn complex simultane- 
ously. This simultaneous bifurcation of an infinite num- 
ber of complex eigenvalues at the phase transition point 



has never been reported before for any PT-symmetric po- 
tentials to our best knowledge 0, ||. Second, we show 
that while stable families of solitons can exist in PT lat- 
tices (below the phase transition point), increasing the 
gain-loss component has an overall destabilizing effect on 
soliton propagation. Specifically, when the gain-loss com- 
ponent increases, the parameter range of stable solitons 
shrinks as new regions of instability appear. Thirdly, we 
investigate the evolution of unstable PT solitons under 
perturbations, and show that the energy of these per- 
turbed solitons can grow unbounded even if the PT lat- 
tice is below the phase transition point. 



II. SIMULTANEOUS COMPLEX-EIGENVALUE 
BIFURCATION AT THE PHASE TRANSITION 
POINT 
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We begin by investigating the bifurcation of the con- 
tinuous spectrum (Bloch bands) of the linear Schrddinger 
operator in Eq. at the phase transition point. 

The phase transition point is a point above which the 
spectrum is no longer purely real as the strength of 
the imaginary (gain-loss) contribution relative to the 
real (refractive-index) contribution in the potential is in- 
creased. We will show that at the phase transition point, 
an infinite number of Bloch bands turn complex simulta- 
neously. 

We first consider this bifurcation in one dimension. In 
this case, the linear Schrddinger equation is 



(2.1) 



iU z + U xx + V(x)U = 0, 



where the PT lattice potential V(x) is given in Eq. (|1.2[) . 
The continuous spectrum of this Schrddinger equation 
consists of Bloch modes of the form 



U(x,z) = p(x;k)e ikx - i » z , 



(2.2) 



where p(x; k) is a 7r-periodic function in x, k is the 
wavenumber in the irreducible Brillouin zone — 1 < k < 
1, and /x is the propagation constant. The values of /i 
and k are related. The relation \i = /i(fc) is called the 
diffraction relation, and all admissible values of /i form 
the continuous spectrum of Eq. (|2.1[) . 

For the PT lattice (|1.2I) . the phase transition point is 
known to be Wo — 0.5 [8]. Below this phase transition 
point (Wo < 0.5), the continuous spectrum is all real 
and comprises an infinite number of segments (known as 
Bloch bands). The gaps between these Bloch bands are 
called bandgaps; the largest, which contains everything 
to the left of the continuous spectrum, is the semi- infinite 
gap and further gaps are numbered (in our case from left 
to right). As an example, at Wo — 0.4 and Vo = 6, 
the diffraction relation is shown in Fig. [TJ and the Bloch 
bands and bandgaps are shown in Fig. [5] 

As Wo increases, bandgaps shrink (see Fig. [5]). At the 
phase transition point (Wo = 0.5), all Bloch bands merge 
(see Figs. [I] and [3]). Above the phase transition point 



FIG. 1: (Color online) Diffraction relations of PT lattices 
(|1.2p for three Wo values 0.4, 0.5 (upper panel) and 0.6 (lower 
panel) at Vo = 6. The inset in the lower right panel is ampli- 
fication of the small boxed region near k — 1 and Im[/i] = 
of the same panel. 
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FIG. 2: (Color online) Bandgap structure of the PT lat- 
tice (|1.2p as Wo crosses the phase transition point 0.5 (with 
Vo = 6). Above this phase transition point, complex eigen- 
values /j, bifurcate out simultaneously from points A, C, ... 
where Bloch bands merge (see the upper panel). The real 
and imaginary parts of these complex eigenvalues versus Wo 
at the Brillouin edge k = 1 are plotted in the upper and lower 
panels respectively. 
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(Wo > 0.5), complex eigenvalues appear in the Bloch 
bands. This phase transition has been reported before 
0. For example, the diffraction relation at Wq = 0.6 
and Vb = 6 is displayed in Fig. Q] It is seen that complex 
eigenvalues /i arise in the Bloch bands near edges k = ±1 
of the Brillouin zone. 

What was not known about this phase transition, how- 
ever, is that right above this phase transition point, 
complex eigenvalues appear simultaneously in an infinite 
number of Bloch bands. To demonstrate, the dependence 
of eigenvalues \i on Wo at Brillouin-zone edge k = 1 is 
shown in Fig. [2j We can see that at the phase transi- 
tion point Wq = 0.5, complex eigenvalues [i bifurcate out 
simultaneously from point A where the first and second 
Bloch bands merge, and from point C where the third 
and fourth Bloch bands merge, with both band-mergings 
occurring at the Brillouin-zone edges k = ±1 (this bifur- 
cation of complex eigenvalues does not occur from point 
B where the second and third Bloch bands merge at the 
Brillouin-zone center k = 0). 

Below, we show analytically that at the phase tran- 
sition point, an infinite number of complex eigenvalues 
bifurcate out simultaneously from an infinite number of 
Bloch bands. In particular, bifurcations of these complex 
eigenvalues occur at points where the (2n — l)-th and the 
2rt-th Bloch bands merge (at k — ±1), but not at points 
where the 2n-th and the (2n + l)-th Bloch bands merge 
(at k = 0), for all positive integers n — 1,2,3,- •■ (see 
Figs, ffland |2}. 

We look for solutions to Eq. (|2.ip of the form U = 
u(x)e~ ltiZ , where u satisfies the equation 



fill 



+ u xx + Vq (cos 2 x + \Wq sin 2x) u = 0. 



(2.3) 



At the phase transition point Wq = 0.5, Eq. (|2.3[) reduces 
to 



M , ,V , o ix . 



f,+ ^-\ u + u xx + ^(e Zlx )u = 0. (2.4) 



Under the variable transformation £ = iWVo/2 e lx , this 
equation becomes Bessel's equation, 

C 2 «« + f«£ + (V -|i-yW 0, (2.5) 



thus it has exact solutions in terms of Bessel functions 




(2.6) 



where k = ± 



V 



+ fc 2 



(2.7) 



This is the exact diffraction relation at the phase tran- 
sition point, as can be seen by utilizing the power-series 



expansion of the Bessel function to expand the above 
Bloch solution (|2.6[) into a Fourier series 




i+k)x 



(2.8) 



where k is seen to be the wavenumber and p(e 2lmx ) is a 
7r-periodic function. By factoring out the 7r-periodic term 
e 2mx f rom e ikx ^£ or a cer t a i n integer n) and combining 
it with p(e 2lmx ), one can restrict the wavenumber k to 
be in the Brillouin zone — 1 < k < 1, as is customary in 
the Bloch theory (see Fig. [T}. The diffraction relation 
(I2.7[) shows that the continuous spectrum at the phase 
transition point Wq — 0.5 is — Vb/2 < fi < oo and is 
entirely real. When k = n is an integer, the two Bessel 
solutions J±k(x) in (|2.6[) are linearly dependent. This 
corresponds to the points where different Bloch bands 
merge (see points A, B, C, ... in Fig. [2}, and the associ- 
ated /x values are 



Vb 



n = 0,1,2, 



(2.9) 



These fi values are located at either fc = 0orfc = ±lof 
the Brillouin zone on the diffraction curves, depending 
on whether n is even or odd (see Fig. [1}, and their Bloch 
functions are 7r-periodic for even n and 27r-periodic for 
odd n. 

We now consider the case where Wo is near the phase 
transition point 0.5, i.e., Vb(Wo — 0.5) = e <C 1. In this 
case, Eq. (|2.3p becomes 



2 



U + Uxx + 



111 

2 



(e 2ix ) u + eisin(2x) = 0, (2.10) 



whose solutions and the corresponding diffraction re- 
lation /i = /i(fc) can be derived by the perturbation 
method. For simplicity, we only derive its solutions u(x) 
which are ir- or 27r-periodic (these Bloch solutions are 
degenerate). The corresponding fi values are then those 
with k = or k = ±1 on the diffraction curves (see Fig. 
[I}. These solutions and the associated fi values can be 
expanded as power series in e 1 / 2 , 



'0.5, 1/J . ,3/2 
M = 7T + n + e n l + en 2 + 6 n 3 



u(x) = u + e 1 / 2 ^! + eu 2 + e' yz ui + 



3/2, 



(2.11a) 
(2.11b) 



where no = 0, 1, 2, • ■ ■ , and coefficients Hi, nz, 113. ■ • • in 
(|2.11a[) are certain constants. Details of this perturbation 
calculation are presented in Appendix 1. The main re- 
sults for these coefficients nj., 112, • • • at various no values 
are summarized in the following table. 

We see from this table that when no = 1, 3, which cor- 
respond to points A, C in Fig. [2J the coefficient ni or 
n3 is imaginary, thus complex eigenvalues bifurcate out 
simultaneously above the phase transition point (e > 0). 
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TABLE I: Coefficients in the fi expansion (12. Hall , 
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This diffraction relation (|2 . 14[) shows that complex eigen- 
values appear in this 2D PT lattice if and only if com- 
plex eigenvalues appear in the ID PT lattice (|1.2p . Thus 
all eigenvalues in the 2D system (|2.12l) are real when 
Wq < 0.5, and a phase transition occurs at Wo = 0.5 
above which complex eigenvalues arise. In addition, an 
infinite number of Bloch bands turn complex simultane- 
ously right above this phase transition point. 

III. STABILITY OF PT SOLITONS IN ONE 
DIMENSION 



In addition, the imaginary part of these complex eigen- 
values at n = 3 is much smaller than that at n = 1 
since the former is of order e 3 / 2 while the latter is of or- 
der e 1 / 2 . However, no complex eigenvalues appear when 
n = 0, 2 (the latter corresponds to the point B in Fig. 
[2]). All these analytical results are in complete qualitative 
and quantitative agreement with Fig. [5] as we have care- 
fully checked. Continuing these calculations to higher no 
values, we have found that the coefficient n,2ro+i is al- 
ways imaginary for no = 2m + 1, where m — 0, 1, 2, • • • . 
Thus complex eigenvalues bifurcate out simultaneously 
from all odd values of n at the phase transition point 
Wo = 0.5. 

The above table also shows that below the phase tran- 
sition point (Wq < 0.5, or e < 0), the eigenvalue /i from 
the expansion (|2.11aj) is always real for all integers no. 
In addition, a gap opens at the corresponding fj, values of 
— Vb/2 + rig. Furthermore, the width of the n th gap is of 
order e 11 / 2 . Above the phase transition point, the even- 
numbered bandgaps reopen, whereas the odd-numbered 
bandgaps close and complex eigenvalues bifurcate out. 
All these analytical conclusions match perfectly with Fig. 
[5] as well. 

Now we consider eigenvalue bifurcations in two dimen- 
sions. In this case, the 2D linear Schrodinger equation 
(PD is 



iU z + U xx + U yy + V(x, y)U = 0, 



(2.12) 



where the PT lattice potential V(x,y) is given in Eq. 
(I1.3[) . This 2D potential is separable, thus the Bloch 
modes of Eq. (|2~T2l are [H[ 

U(x, y, z) = e ^+^y-^p( x; kl )p(y; k 2 ), (2.13) 

where p(x; k) is the ID 7r-periodic function as given in 



/i = ju(fci) + p,(k 2 ) 



(2.14) 



is the 2D diffraction relation, k\ , k 2 are Bloch wavenum- 
bers in x and y directions and are located inside the irre- 
ducible Brillouin zone — 1 < k%,k2 < 1, and the function 
ju(fc) is the diffraction relation of the ID equation (|2.1I) . 



In the presence of cubic nonlinearity, the mathematical 
model becomes the NLS equation (jl.lj) with a PT lattice 
potential. In this case, light can self-localize and form 
solitons. In this section, we study these PT solitons and 
their linear-stability behaviors in one dimension. 

In one dimension, the NLS equation (jl.ip becomes 



iU z + U xx + V(x)U + a\U\ 2 U = 0. 



(3-1) 



Here the PT lattice V(x) is taken as (TPl) with Vo = 6, 
and a = dbl. Solitons in this model are sought of the 
form 



U(x,z)=e- l »*u(x), 



(3.2) 



where u(x) is a localized function, and /i is a real propaga- 
tion constant. These solitons can be computed by either 
the squared operator iteration method or the Newton- 
conjugate-gradient method applied to the normal equa- 
tion [15| . They exist when /i lies inside bandgaps of 
the linear system for Wo both below and above the 
phase transition point. Above the phase transition point 
(Wo > 0.5), linear waves amplify exponentially during 
propagation, thus any solitons would also be unstable to 
perturbations. So we only need to consider Wo < 0.5 
below. 

To determine the linear stability of these PT solitons, 
we perturb them as 



U 



ltiz u{x) + u{x)e Xz + w*(x)e 



,A*2 



(3.3) 



where \w\ <c \u\. After substitution into equation 
(|3.1|) and linearizing, we arrive at the eigenvalue problem 





where 



£ = 



= A 



\ L 2 i L22 J 
Ln=M + d xx + V(x) + 2a\u\ 2 1 

L\2 = (TV 2 , 

L 2 i = -a (u 2 )* , 
L22 = 



(3.4) 



(v + d xx + V*{x) + 2a\u\ 2 ) 



5 



This eigenvalue problem can be computed by the Fourier 
collocation method (for the full spectrum) or the Newton- 
conjugate-gradient method (for individual discrete eigen- 
values) [15(- If eigenvalues with positive real parts exist, 
the soliton is linearly unstable; otherwise it is linearly 
stable. 

We first consider PT solitons in the semi-infinite gap 
under focusing nonlinearity (a — 1). For Wq = 0.45, 
two families of PT solitons are obtained and their power 
curves are displayed in Fig. [3] (left). Here the power of a 
soliton is defined as 



10 



J — ! 



\u(x] fi)\ 2 dx. 



(3.5) 



In this figure, the lower power curve is for the funda- 
mental solitons which exhibit the same PT symmetry 
u* (x) = u(—x) and whose real parts possess a single dom- 
inant peak. The profile of such a soliton at [i — —3.5 is 
displayed in Fig. [3] (right). This soliton family bifurcates 
out of the first Bloch band, and the solitons near this 
Bloch band are low-amplitude Bloch-wave packets. We 
have found that the entire branch of this fundamental- 
soliton family is linearly stable, which is indicated by 
solid lines of its power curve in Fig. [3] (left). The up- 
per power curve in Fig. [3] is for the dipole solitons. 
This power curve features double branches which ter- 
minate before reaching the first Bloch band (a similar 
phenomenon occurs in purely real lattices [15l. Il6|). Pro- 
files of three such solitons on the lower power branch 
are displayed in Fig. |4] (top). It is seen that the real 
parts of these dipole solitons possess two dominant peaks 
of opposite phase (which is why they are termed dipole 
solitons). Unlike the fundamental solitons, these dipole 
solitons are linearly stable only in a certain portion of 
their existence region. Specifically, only dipole solitons 
on the lower branch and with /j, < —3.8 are stable (see 
Fig. |3j (left)). For dipole solitons in this region, their 
spectra are entirely imaginary (see Fig. |4] (bottom left)). 
At /i = —3.8, stability switching occurs where a quadru- 
ple of complex eigenvalues bifurcate off of the edge of 
the continuous spectrum (see Fig. |4] (bottom center)). 
Within the unstable region, there is a second eigenvalue 
bifurcation at fj, » —3.4 of the lower branch (near and 
on the left side of the power minimum) where a pair of 
real eigenvalues bifurcate from zero (see Fig. |4] (bot- 
tom right)). Some of these stability behaviors on dipole 
solitons are similar to those in the purely real potential 
(Wo = 0) [Hj]. A notable difference is that for real poten- 
tials real eigenvalues bifurcate out of the origin exactly 
at the minimum of the power curve [15|, whereas here 
this real-eigenvalue bifurcation occurs not at the power 
minimum. An analytical explanation for this new phe- 
nomenon will be given in Appendix 2. 

Next we consider PT solitons in the first gap under 
defocusing nonlinearity (er = —1). Again, for Wo = 0.45, 
two families of PT solitons are obtained and their power 
curves are displayed in Fig. [5] (left) with stability re- 
sults indicated. The lower curve is for fundamental soli- 
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FIG. 3: (Color online) One-dimensional PT solitons in the 
semi-infinite gap under focusing nonlinearity (a — 1) for Vo = 
6 and Wo = 0.45. (left) Power curves of these solitons; the 
lower curve is for fundamental solitons and the upper curve 
for dipole solitons; solid and dashed lines represent stable and 
unstable solitons respectively (the same holds for all other 
figures); the shaded region is the first Bloch band, (right) 
Profile u(x) of a fundamental soliton at fi = —3.5 (marked 
by a dot on the lower curve of the left panel); the solid blue 
line is for the real part and dashed pink line for the imaginary 
part. 
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FIG. 4: (Color online) Dipole solitons (top) and their linear- 
stability spectra (bottom) for three /x values in the semi- 
infinite gap. The power curve of these dipole solitons is shown 
in Fig. [3] and the locations of these solitons are marked by 
dots on that power curve. 



tons whose profiles at two /z values are depicted in Fig. 
[5] (top right panel), while the upper curve is for dipole 
solitons, whose profiles are similar to those in Fig. [5] 
(middle panel) below. The fundamental-soliton family 
bifurcates out of the first Bloch band, whereas the dipole 
family docs not. We have found that all solitons in this 
dipole family are linearly unstable (see Fig. [5] (left)). 
The fundamental-soliton family, however, is linearly sta- 
ble when [x < —1.77. At /i = —1.77, stability switching 
occurs where a pair of real eigenvalues bifurcate out from 
zero (see Fig. [3]). Notice that unlike in real potentials 
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[15| . this zero-eigenvalue bifurcation does not occur at a 
power extremum since the potential here is complex. An 
explanation for this will be presented in Appendix 2. 



15 



is caused by a quadruple of complex eigenvalues (see Fig. 
[5] (bottom panel)). 
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FIG. 5: (Color online) One-dimensional PT solitons in the 
first gap under defocusing nonlinearity (a = — 1) for Vo = 6 
and Wo = 0.45. (left) Power curves of these solitons; the 
lower curve is for fundamental solitons and the upper curve 
for dipole solitons; (top right) two fundamental solitons at 
fi = —2 and —1.7 (marked by dots in the left panel); (bottom 
right) linear-stability spectra of these solitons. 

The stability results of PT solitons in Figs. [3]to[S]were 
obtained for a specific Wo value of 0.45. Now we discuss 
how these stability results change when Wq steadily in- 
creases from to 0.5. First we consider PT solitons in 
the semi-infinite gap under focusing nonlinearity. We find 
that when < Wq < 0.5, the entire fundamental-soliton 
family remains stable. The dipole family, however, is 
stable only on the left side of its lower branch (see Fig. 
02), and this stable region shrinks as Wq increases. Next 
we consider PT solitons in the first gap under defocusing 
nonlinearity. When Wq = (i.e., the lattice is real), the 
fundamental-soliton family is stable in the entire first gap 
[l5j |. As Wo rises above 0.3, an unstable region grows off 
the edge of the second band. At the stability switching 
point a pair of real eigenvalues bifurcate from zero as il- 
lustrated in Fig. [SJ Regarding the dipole-soliton family, 
its entire upper branch is unstable for all Wq values. Its 
entire lower branch is also unstable when Wo > 0.44. For 
Wo < 0.44, a certain portion of its lower branch is stable; 
but as Wo increases, this stable region shrinks and then 
totally disappears when Wo > 0.44. To demonstrate this 
reduced stability of dipole solitons with increasing Wo, 
the power curves of these dipole solitons at two Wo val- 
ues of 0.35 and 0.4 are shown in Fig. [5] (top panel) with 
stability results indicated. The soliton profiles at /i = — 2 
are also shown in the middle panel of the same figure. It is 
seen that the stable region of dipole solitons at Wo =0.4 
is much shorter than that at Wo = 0.35. Notice also that 
as Wo increases the width of the first gap decreases which 
is often a sign of decreased stability. The unstable region 
on the lower branch is largely located near the edge of 
the second Bloch band, and the instability in this region 
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FIG. 6: (Color online) Dipole solitons in the first gap under 
defocusing nonlinearity (cr = — 1) at Wo values 0.35 and 0.4 
(with Vo = 6). (top) Power curves; (middle) soliton profiles 
at fi — —2 (marked by dots in the top panel); (bottom) linear- 
stability spectra of the solitons in the middle panel. 

The above stability results of PT solitons show that as 
Wq increases (but still below the phase transition point), 
the stable regions of PT solitons generally shrink (see Fig. 
[5]). The only exception is the fundamental-soliton fam- 
ily in the semi-infinite gap under focusing nonlinearity, 
which remains entirely stable up to the phase transition 
point. Overall, the inclusion of the gain-loss component 
in the PT lattice has a destabilizing effect on solitons. 



IV. STABILITY OF PT SOLITONS IN TWO 
DIMENSIONS 

In this section we analyze the linear stability of solitons 
in a two-dimensional PT lattice. We will show that the 
destabilizing effect of the gain-loss component is more 
prominent in this case, even for fundamental solitons in 
the semi-infinite gap. 

In two dimensions, the mathematical model is Eq. 
(ft~T|) . or 

iU z + U xx + U yy + V(x, y)U + <r\U\ 2 U = 0, (4.1) 
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where the PT lattice V (x, y) is taken as (| 1 .3[) with Vq = 6. 
Solitons in this model are sought of the form 

U(x,y,z)=c- i ^u{x, V ) 1 (4.2) 

where u(x, y) is a localized function, and /i is a real prop- 
agation constant. These solitons as well as their linear- 
stability spectra can be obtained by numerical methods 
similar to the ID case. The phase transition point in 
this 2D model is also Wo = 0.5, above which all solitons 
are linearly unstable. Thus we only consider Wq < 0.5 
below. 

For simplicity we only consider 2D fundamental PT 
solitons in the semi-infinite gap under focusing nonlin- 
earity (er = 1). These fundamental solitons possess the 
PT symmetry u* (x, y) — u(—x, —y), and their real parts 
have a single dominant peak. Profiles of such solitons can 
be found in Fig. [8] (upper panel) later. We find that these 
fundamental solitons are stable only on a finite /i-interval 
even for small values of Wq. In addition, this stable re- 
gion shrinks as Wo increases and totally disappears when 
Wo > 0.47. To demonstrate, power curves of these soli- 
tons as well as their stability regions at two Wo values 
of 0.2 and 0.3 are displayed in Fig. [7J It is seen that 
the stable region is finite even though the existence re- 
gion of solitons is infinite. In addition, as Wo increases 
from 0.2 to 0.3, the stable region has shortened by sev- 
eral times. For each Wq, there are two unstable regions, 
one located at large negative fi values, and the other one 
located near the first Bloch band. For large negative val- 
ues of fx the instability is due to a quadruple of complex 
eigenvalues, whereas for /i values near the first band, the 
instability is due to a pair of real eigenvalues. Examples 
of the spectrum in each region are shown in Fig. [5] with 
Wo = 0.3. We see that in this 2D case, even the funda- 
mental solitons in the semi-infinite gap are destabilized 
by the addition of the gain-loss component in the lattice. 



W Q = 0.3 



fj, = -10 



H = -8.5 



H = -7.3 





FIG. 7: (Color online) Power curves of fundamental 2D 
solitons in the semi-infinite gap under focusing nonlinearity 
(a = 1) for Vo = 6 and two Wo values of 0.2 and 0.3. The 
inset in the right panel is amplification of the power curve 
near the first Bloch band in the same panel. 
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FIG. 8: (Color online) Fundamental 2D solitons (|u(x, 
(top) and their linear-stability spectra (bottom) for three fi 
values in the semi- infinite gap with a — 1, Vo = 6 and Wo = 
0.3. The power curve of these solitons is shown in Fig. 
(right panel). 



V. NONLINEAR EVOLUTION OF PT 
SOLITONS UNDER PERTURBATIONS 



In this section, we examine the nonlinear evolution 
of PT solitons under weak perturbations. We find that 
when a PT soliton is linearly stable, then it is also non- 
linearly stable and propagates robustly against perturba- 
tions. If the soliton is linearly unstable, then it breaks up 
under perturbations, and its amplitude and energy can 
grow unbounded over distance. 

First we consider the ID fundamental soliton shown 
in Fig. |3l which resides in the semi-infinite gap under 
focusing nonlinearity and is linearly stable. We perturb 
it by 5% random noise perturbations and then simulate 
its evolution in Eq. (|3.1[) . The simulation result is shown 
in Fig. [9] (left). We can see that even after z = 100 units 
of propagation, this soliton remains robust and does not 
break up. Thus this soliton is also nonlinearly stable. 
Next we consider the ID fundamental soliton shown in 
Fig. [5l which resides in the first gap under defocusing 
nonlinearity and is linearly unstable. When this soliton 
is perturbed by 5% random noise perturbations, its evolu- 
tion is shown in Fig. |H] (right). It is seen that this soliton 
quickly blows up and spreads out, thus is obviously non- 
linearly unstable. Notice that the peak amplitude and 
energy of this perturbed soliton steadily increase with- 
out bound over distance. This indicates that the gain- 
loss component of the PT lattice steadily feeds energy 
into the solution. Recall that the Wo value in this case is 
below the phase transition point, thus linear waves do not 
grow. Consequently the energy growth in this evolution 
is solely due to the nonlinear effects. 

Lastly we consider the 2D fundamental soliton shown 
in Fig. [5] (left panel) , which resides in the semi- infinite 
gap under focusing nonlinearity and is linearly unstable. 
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a = 1 fi = —3.5 a = — 1 fx = —1.7 




FIG. 9: (Color online) (left) Nonlinear evolution of the stable 
ID soliton in Fig. [3] under 5% random noise perturbations; 
(right) Nonlinear evolution of the unstable ID soliton in Fig. 
[5] (with n = —1.7) under 5% random noise perturbations. 
Shown is the field \U(x, z)\ in the (x, z) plane. 



increasing the gain-loss component of the lattice has an 
overall destabilizing effect on soliton propagation. In ad- 
dition, we have shown that when unstable PT solitons 
are perturbed, the energy of the solution can grow un- 
bounded even though the PT lattice is below the phase 
transition point. 
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Appendix 1: Calculation of eigenvalue bifurcations 
at the phase transition point 



When this soliton is perturbed by 5% random noise per- 
turbations, its evolution is shown in Fig. 1101 It is seen 
that the power (and peak amplitude) of this perturbed 
soliton also grows oscillatorily without bound, thus this 
soliton is nonlinearly unstable. This oscillatory growth 
occurs since the unstable eigenvalues of this soliton are 
complex (see Fig. [8] (lower left panel)). 
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FIG. 10: (Color online) Nonlinear evolution of the unsta- 
ble 2D soliton in Fig. [8] (with /i = —10) under 5% random 
noise perturbations, (left) Power evolution versus distance z; 
(right) Solution profiles at two distances z = 6.6 and 7.3. 



VI. SUMMARY 

In summary, we have analyzed the linear phase tran- 
sition and nonlinear solitons in PT-symmetric photonic 
lattices. We have shown that at the phase transition 
point, an infinite number of linear Bloch bands turn com- 
plex simultaneously. We have also shown that while con- 
tinuous ranges of stable solitons can exist in PT lattices, 



In this appendix, we calculate eigenvalue bifurcations 
at the phase transition point in Eq. (|2 . 10[) by pertur- 
bation methods. The solution u(x) to this equation is 
required to be it- or 27r-periodic, and perturbation ex- 
pansions for u(x) and eigenvalue /i are as given in Eq. 

Let us define the operator 



L — d x 



V 



2ix 



(A.l) 



After substituting expansions (|2.11[) into Eq. (|2.10[) and 
collecting terms of the same order in e 1 ' 2 we arrive at the 
following system of linear equations 



Lu 
Lui 



0. 



-niu , 

-isin(2a;)u m _2 



n-jUm— j 

3 = 1 



(A.2a) 
(A.2b) 

(A.2c) 



for m — 2, 3, 4, • • • . The solution u is 



where 



(Yo/sy 



o. rn c 



i(2m+no)x 



m!(m + no)! 



for m > 0, 



(A.3) 



(A.4) 



and a m = for m < 0. This solution comes directly 
from (|2.8I) by replacing the wavenumber k with the inte- 
ger no. The remaining linear inhomogeneous equations 
()A.2b[) - (|A.2c[) for ui,u 2 ,.-. will be solved by first im- 
posing the solvability condition due to the Fredholm Al- 
ternative Theorem and then expanding the solution in 
terms of Fourier series. 
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The adjoint operator of L is 



(A.5) 



and has kernel Uq since L a Uq = 0. The Frcdholm Alter- 
native Theorem requires that the forcing terms in equa- 
tions (|A.2I) be orthogonal to Uq. As we have mentioned 
earlier, we are concerned with n- and 27r-periodic solu- 
tions here and thus define the inner product as 

(f(x),g(x)) = ±- [' f(x)g*(x)dx. (A.6) 

27r J-7T 

Using the fact that (e ipx , e lqx ) — S Piq for integers p, q, we 
obtain the inner product 



When no = (hence ni = u\ = 0), its solvability condi- 
tion is 



n 2 



(isin(2a;)uo, Uq) 
(«o,«$) 



(A.14) 



which gives n 2 = Vq/8. For no > 1, after substituting in 
the solution (IA.9I) - (|A.10I) for u±, the solvability condition 
of (|A.13|) gives 



(i sui(2x)uq , Uq) 

(L _1 MO,Uo) 



(A.15) 



By rewriting isin(2x) = | (e l2x — e~ l2x ) we may again 
use the orthogonality of the Fourier modes to work out 
the inner products explicitly, 



0, n = 1,2,3, 



(A.7) 



Eq. (IA.2b[) for u\ has the solvability condition 

= -m(«o,«S>. (A.8) 



(i sin(2a;)wo, Wq) = — -Oq, for no = 1, (A. 16a) 



(i sin(2a;)wo, Wq) = 0, 



for n =2,3,4..., 

(A.16b) 



(L 1 u ,Uq) = &- no ao, for n = 1,2,3.. 



(A.16c) 



Thus, when no = then ni = 0. For other no this 
solvability condition is satisfied automatically and the 
solution u\ may be formally written as 



Thus, 



Ml = — ni-L Uq. 

Expanding L _1 mq into Fourier series 



L 1 u a = ^2 bme 1 



(2m+no)s 



(A.9) 



(A.10) 



m— — oo 



and substituting it into L \L 1 uo] = uo we find that the 
coefficients b m satisfy the recursion relation 

- 4 (m 2 + mn ) b m + -y&m-i = a m (A. 11) 

for all integers m. The relevant coefficients are 



b-i = —a 
Vo 

b-2 = -T72"( n o - l)ao 



(A.12a) 
(A.12b) 



(_iyio-l / 8 \ n ° 

6-n = ^4 (n ° _1)!2 (^J 00 (AJ2C) 



b m = 0, for m < — n . 



(A.12d) 



Notice that this series also terminates in the negative m 
direction at m = —no. 

The equation (IA.2c|) for u 2 is 



Lu2 = — isin(2a;)uo — niui — n 2 UQ. 



(A.13) 



K 1/2 

ni = ±i — - — , for no = 1, 



(A.17) 



and ni = for no > 1. This means if no = 1 then ni is an 
imaginary number and, returning to the expansion for \x 
in equation (|2.11al) . that fi is a complex number for Wq 
above the phase transition point, e > 0, and real below, 
e < 0. This is the bifurcation that occurs at edge of the 
Brillouin zone where the first and second bands merge 

(see Figs. [Hand O. 

For n = 1, we can proceed to solve Eq. (|A.13j) for m 2 
by Fourier expansion. Then from the solvability condi- 
tion for the U3 equation we can find that 

ii2 = , for n = 1. 
For no > 2 we formally write the solution u 2 as 

u 2 = —L^ 1 [i sin(2a;)uo] — n 2 £ t«o (A. 18) 

since we know that L _1 [i sin(2a;)uo] is well defined in view 
of the orthogonality (|A.16b|) . Expanding it into a Fourier 
series 



L 1 [i sin(2a:)uo] = } c m e 



m— — oo 



(A.19) 



it is easy to find that the coefficients c m satisfy the re- 
cursion relation 



V 



1 



4 (ra 2 + mn ) c m + —Cm-i = 77 - a m +i) , 



(A.20) 
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and 



C-l 



1 2 

-ao, C-2 = -7— 7TT a O- 



(n + 1) 



(1 + n )Vo 



Again there are only a finite number of terms in the neg- 
ative m direction, i.e. c m = for m < —no. 

For no > 2 (hence ni = U\ = 0), Eq. (|A.2cj) for 1*3 is 



I/M3 = -n 3 ii , 



thus 



(A.21) 
(A.22) 



U3 = -n 3 L 1 u Q , 
and Eq. (IA.2cp for U4 is 

L114 — — isin(2x)u 2 — n 2 u 2 — n 4 ?io- (A. 23) 



After substituting in (|A.18I) for m 2 the solvability condi- 
tion is 

= nl(L^ 1 u Q ,u^ t ) 

+ 2n2(L _1 isin(2a;)wo, Uq) 

+ (isin(2x)L _1 isin(2a;)uo,M5}, (A. 24) 

with coefficients given by 

{L~ x u ,ul) = &- no ao, (A.25a) 
(L _1 isin(2x)u ,-u5) = c_ no a , (A. 25b) 

(isin(2x)L _1 isin(2x)uo,Wo) = -- (c_ no ai + c_ no+1 a ) . 

(A.25c) 

For n = 2 this gives n 2 = V a /2A, -5Vb/24. For n > 3 
we find that n 2 is a double root, 



V Q 1 

n 2 = - — -2~ 



8^-1 



(A.26) 



At these n 2 values, the solution U4 is well defined and is 
given by 

U4 = — [isin(2x)M 2 + n 2 w 2 ] — rx^L - uq. (A. 27) 
For n > 2, Eq. (|A~2cl) for «5 is 
Lus = — isin(2a;)M3 — n-3it 2 — n 2 it3 — nsitQ. (A. 28) 



After substituting in equations (|A.18[) and (| A.22I) the 
solvability condition for this W5 equation reduces down 
to 

= n 3 [n^L-^oX) + (i _1 isin(2x)u , u*)] . (A.29) 

Thus, for no = 2 we must have n3 = 0; and for no > 3 
this condition is satisfied automatically since n 2 from Eq. 
(|A~26t is a double root o f Eq. (|A34| . 
For n > 3, Eq. (|A~2cl for Uq is 

Lue = — i sin(2a;)u4 — n 2 U4 — n3^3 — n4M 2 — ngito. (A. 30) 



Substituting in Eq. (|A.22I) and noting that (m 2 ,Uq) = 
(in view of (|A.29p ) we are left with the solvability condi- 
tion 



(«4(n 2 + i sin(2a;)), Uq) 



6 (L-%o,^) 
This condition may be further simplified, 

(u 4 (n 2 +isin(2a;)),Mo) 

= (u 4 , (n 2 - isin(2x))Mg) 

= -(isin(2a;)u2 + n 2 u 2 - n 4 u , u 2 ) 

= — (i sin(2x)u 2 , w 2 ). 



(A.31) 



Thus for no = 3, we get 



113 



±i y ° 



3/2 



2 9 



(A.32) 



and for no > 3, we get n3 = 0. This shows that there 
is another bifurcation point of complex eigenvalues at 
n = 3, where the third bandgap closes (see Figs. [T]and 

The results of the above perturbation calculations are 
summarized in Table 1 of the main text. Continuing 
these calculations to higher no values, we have found that 
the coefficient n 2m+ i is always imaginary for no = 2m+l, 
where m = 0, 1, 2, • • • . Thus complex eigenvalues bifur- 
cate out simultaneously from no = 1, 3, 5, • • • at the phase 
transition point Wq = 0.5. 



Appendix 2: Analytical criterion for zero-eigenvalue 
bifurcation of solitons in complex potentials 

In real potentials (such as when Wo — in fll.2[) ). the 
power curve does more than just a convenient way to cat- 
alogue and parameterize a continuous family of solitons 
for various values of the propagation constant fi. Specifi- 
cally, whenever the power curve has a local extremum the 
zero eigenvalue in the linear-stability spectrum of solitons 
then bifurcates out along the real and imaginary axes on 
the two sides of the power extremum respectively, lead- 
ing to a change of stability at the pow er extremum (if no 
other unstable eigenvalues exist) jla |. In this appendix 
we consider the extension of this concept to general com- 
plex potentials (which include PT-symmetric lattices as 
special cases) . The resulting analytical criterion for zero- 
eigenvalue bifurcation will explain the stability switch- 
ings in Fig. [5] and Fig. [7] (right side), as well as the onset 
of real eigenvalues in Fig. 2] (right panel). 

Let us begin with the eigenvalue problem (|3.4p derived 
in the main text, 




(A.33) 



where we know that A = is always an eigenvalue with 
algebraic multiplicity of at least two due to phase invari- 
ance of Eq. (|l.ll) . The eigenfunction and generalized 
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eigenfunction of this zero eigenvalue associated with the 
phase invariance can be written explicitly in terms of the 
soliton u(x), 



Here the superscript "T" stands for transpose of a ma- 
trix. Then the solvability condition of Eq. (|A.34j) is 



U 

-u* 



and C 



u 

-u* 



^ = 0, 



(A.37) 



Thus, for nonzero eigenvalues to bifurcate out from the 
origin, A = must have algebraic multiplicity of at least 
3 at that point. A sufficient condition for this to occur is 
that there be a second generalized eigenfunction -0 which 
solves 



which is a sufficient condition (criterion) for zero- 
eigenvalue bifurcation in general complex potentials. 

For real potentials, it is easy to see that 



(A.34) 



We now use the Fredholm Alternative Theorem to derive 
the solvability condition for Eq. (|A.34[) . Let us denote 
the kernel of the adjoint operator C A as 4>^ A \ i- e -j 



C A 4>^ = 0, 
where the adjoint operator is 



(A.35) 



(A.36) 




(A.38) 



thus the above criterion reduces to P'(n) = 0, i.e., the 
extremum of the power curve [l5j . For general complex 
potentials, however, (f>( A ^ is not equal to the above ex- 
pression, thus stability switching will no longer occur at 
a power extremum. An example of this has been seen in 
Fig. [3 
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